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ABSTRACT 

Blazars have been suggested as possible neutrino sources long before the re¬ 
cent IceCube discovery of high-energy neutrinos. We re-examine this possibil¬ 
ity within a new framework built upon the blazar simplified view and a self- 
consistent modelling of neutrino emission from individual sources. The former 
is a recently proposed paradigm that explains the diverse statistical properties 
of blazars adopting minimal assumptions on blazars’ physical and geometrical 
properties. This view, tested through detailed Monte Carlo simulations, repro¬ 
duces the main features of radio, X-ray, and 7 -ray blazar surveys and also the 
extragalactic 7 -ray background at energies > 10 GeV. Here we add a hadronic 
component for neutrino production and estimate the neutrino emission from BL 
Lacs as a class, “calibrated” by fitting the spectral energy distributions of a pre¬ 
selected sample of BL Lac objects and their (putative) neutrino spectra. Unlike 
all previous papers on this topic, the neutrino background is then derived by 
summing up at a given energy the fluxes of each BL Lac in the simulation, all 
characterised by their own redsliift, synchrotron peak energy, 7 -ray flux, etc. 
Our main result is that BL Lacs as a class can explain the neutrino background 
seen by IceCube above ~ 0.5 PeV while they only contribute ~ 10% at lower 
energies, leaving room to some other population(s)/physical mechanism. How¬ 
ever, one cannot also exclude the possibility that individual BL Lacs still make a 
contribution at the ~ 20% level to the IceCube low-energy events. Our scenario 
makes specific predictions testable in the next few years. 

Key words: neutrinos — radiation mechanisms: non-thermal - BL Lacertae 
objects: general - gamma-rays: galaxies 


1 INTRODUCTION 

Blazars are a class of Active Galactic Nuclei (AGN), 
which host a jet oriented at a small angle with respect 
to the line of sight. Highly relativistic particles moving 
within the jet and in a magnetic field emit non-thermal 
radiation (Blandford & Rees 1978; Urry & Padovani 
1995). This is at variance with most other AGN whose 
energy is mainly thermal and produced through accretion 
of matter onto a supermassive black hole. Because of their 
peculiar orientation and highly relativistic state, blazars 
are characterised by distinctive and extreme observa¬ 
tional properties, including superluminal motion, large 


* E-mail: ppadovan@eso.org 
f Einstein Postdoctoral Fellow 


and rapid variability, and strong emission over the en¬ 
tire electromagnetic spectrum. The two main blazar sub¬ 
classes, namely BL Lacertae objects (BL Lacs) and flat- 
spectrum radio quasars (FSRQ), differ mostly in their 
optical spectra, with the latter displaying strong, broad 
emission lines and the former instead being characterised 
by optical spectra showing at most weak emission lines, 
sometimes exhibiting absorption features, and in many 
cases being completely featureless. 

The spectral energy distributions (SEDs) of blazars 
are composed of two broad humps, a low-energy and 
a high-energy one. The peak of the low-energy hump 
Epeak) can occur at widely different frequencies, ranging 
from about ~ 10 12 ' 5 Hz (~ 0.01 eV) to ~ 10 18 ' 5 Hz (~ 13 
keV). The high-energy hump, which may extend up to ~ 
10 TeV, has a peak energy that ranges between ~ IO 20 Hz 
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(~ 0.4 MeV) to ~ 10 26 Hz (~ 0.4 TeV)(Giommi et al. 
2012 b; Arsioli et al. 2015). Based on the rest-frame value 
of zZpeak, BL Lacs can be further divided into Low energy 
peaked (LBL) sources (^p ea k < 10 14 Hz [< 0.4 eV]), In¬ 
termediate (10 14 Hz< t'peak < 10 15 Hz [0.4 eV < !Zp ea k< 4 
eV)] and High (t/p ea k > 10 15 Hz [> 4 eV]) energy peaked 
(IBL and HBL) sources respectively (Padovani & Giommi 
1995). It is generally accepted that the low-energy blazar 
emission is the result of electron synchrotron radiation, 
with the peak frequency reflecting the maximum energy 
at which electrons can be accelerated (e.g. Giommi et 
al. 2012b). On the other hand, the origin of their high- 
energy emission is still under debate. In the conventional 
leptonic scenarios, 7 -ray emission is thought to be due 
to inverse Compton radiation (e.g. Maraschi et al. 1992; 
Sikora et al. 1994) whereas in leptohadronic scenarios it 
may be the result of proton synchrotron radiation (Aha- 
ronian 2000; Miicke & Protheroe 2001) or may have a 
photohadronic origin (e.g. Petropoulou et al. 2015, and 
references therein). Given that the presence of relativistic 
electrons is necessary to explain at least the low-energy 
hump of the SED, it is reasonable to assume that pro¬ 
ton acceleration takes also place in blazar jets (Bier- 
mann & Strittmatter 1987; Sironi et al. 2013; Globus et 
al. 2014 and references, therein). Although leptohadronic 
models bear several attractive features from the theo¬ 
retical point of view (e.g. Halzen & Zas 1997), these 
are not sufficient to rule out leptonic models. The most 
promising way to settle the issue of proton acceleration 
in blazars is through high-energy neutrino experiments, 
since neutrino emission is an unavoidable outcome of lep¬ 
tohadronic models. 

Padovani & Resconi (2014) have recently extended 
the domain over which blazars could be relevant astro- 
physical sources into neutrino territory. Namely, on the 
basis of a joint positional and energetic diagnostic, they 
have suggested a possible association between eight BL 
Lacs (all HBL) and seven neutrino events reported by 
the IceCube collaboration (IceCube Collaboration 2014). 
Following up on this idea, Petropoulou et al. (2015) have 
modelled the SEDs of six of these BL Lacs using a one- 
zone leptohadronic model and mostly nearly simultane¬ 
ous data. The neutrino flux for each BL Lac was self- 
consistently calculated, using photon and proton distri¬ 
butions specifically derived for every individual source. 
The SEDs of the sources, although different in shape and 
flux, were all well fitted by the model using reasonable pa¬ 
rameter values. Moreover, the model-predicted neutrino 
flux and energy for these sources were of the same or¬ 
der of magnitude as those of the IceCube neutrinos. In 
two cases, i.e. MKN 421 and H 1914—194, a suggestively 
good agreement between the model prediction and the 
neutrino flux was found. 

The hypothesis put forward by Padovani & Resconi 
(2014) and Petropoulou et al. (2015), if correct, should 
materialise in an IceCube detection but this has not hap¬ 
pened yet. At present, in fact, IceCube has not identi¬ 
fied any point sources and therefore its signal remains 
unresolved, although the published upper limits are still 
not ruling out the scenario described above (Padovani & 
Resconi 2014). 

In this paper we aim to calculate the cumulative neu¬ 


trino emission from all BL Lacs within the leptohadronic 
scenario for their 7 -ray emission in order to see if BL Lacs 
as a class can indeed explain the IceCube detections. 
By considering only the neutrino emission produced in 
the blazar jet, we calculate the neutrino background 1 
(NBG) from blazars. This is the most conservative es¬ 
timate, since we do not take into account the ‘cosmo- 
genic’ neutrino production due to the propagation of es¬ 
caping cosmic rays (CRs) in the interagalactic medium 
(Alders & Halzen 2012). As discussed below, calculating 
the NBG from blazars requires a very detailed knowledge 
of the blazar population in terms of t'peaki 7 -ray fluxes, 
redshift, etc. All these parameters, and many more, are 
available in the simulations done in a series of papers by 
Giommi, Padovani, and collaborators. 

Giommi et al. (2012a) (hereafter Paper I) have pro¬ 
posed a new blazar paradigm where blazars are clas¬ 
sified into FSRQ and BL Lacs according to a vary¬ 
ing mix of the Doppler-boosted radiation from the jet, 
the emission from the accretion disc, the broad-line re¬ 
gion, and the light from the host galaxy. This is also 
based on minimal assumptions on the physical proper¬ 
ties of the non-thermal jet emission, and unified schemes. 
These posit that BL Lacs and FSRQ are simply low- 
excitation (LERGs)/Fanaroff-Riley (FR) I and high- 
excitation (HERGs)/FR II radio galaxies with their jets 
forming a small angle with respect to the line of sight 
(e.g. Urry & Padovani 1995). They called this new ap¬ 
proach the blazar simplified view (BSV). By means of 
detailed Monte Carlo simulations, Paper I showed that 
the BSV scenario is consistent with the complex observa¬ 
tional properties of blazars as we know them from all the 
surveys carried out so far in the radio and X-ray bands, 
solving at the same time a number of long-standing is¬ 
sues. 

In a subsequent paper (Giommi, Padovani, & Po¬ 
lenta 2013, hereafter Paper II) the Monte Carlo simula¬ 
tions were extended to the 7 -ray band (100 MeV - 100 
GeV) and found to match very well the observational 
properties of blazars in the Fermi -LAT 2-yr source cat¬ 
alogue (Nolan et al. 2012; Ackermann et al. 2011) and 
the Fermi -LAT data of a sample of radio selected blazars 
(Giommi et al. 2012b; Planck Collaboration 2011). 

Padovani & Giommi (2015) (hereafter Paper III) 
considered the case of very high energy (VHE) emission 
(E > 100 GeV) extrapolating from the expectations for 
the GeV band, and made detailed predictions for current 
and future Cherenkov facilities, including the Cherenkov 
Telescope Array. Finally, Giommi & Padovani (2015) 
(hereafter Paper IV) have extended the predictions below 
the sensitivity of current surveys and estimated the con¬ 
tribution of blazars to the X-ray and 7 -ray extragalactic 
backgrounds. They found that the integrated light from 
blazars can explain the cosmic background at energies > 


1 We use “background” in the astronomical sense of total 
emission from a population of astrophysical objects. This in¬ 
cludes resolved and unresolved sources, the contribution of 
the latter corresponding to the “diffuse” intensity generally 
referred to by high-energy physicists. 
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10 GeV, and contribute « 40 — 70% of the 7 -ray diffuse 
radiation in the 0.1 — 10 GeV band. 

The purpose of this paper is to merge the simulation 
and theoretical approaches by adding to the former a 
hadronic “prior” necessary for the neutrino emission, as 
detailed in Sect. 3. As in papers I - IV we use a ACDM 
cosmology with Ho = 70 km s -1 Mpc^ 1 , = 0.27 and 
11 a = 0.73 (Komatsu et al. 2011). 


2 THEORETICAL MODELLING 

Neutrino production is studied within a specific theoret¬ 
ical framework for blazar emission where photohadronic 
interactions have an active role in shaping the blazar 
SED, as detailed in Petropoulou et al. (2015). In their 
model, the low-energy emission of the blazar SED is at¬ 
tributed to synchrotron radiation of relativistic electrons, 
whereas the observed high-energy (GeV - TeV) emission 
has a photohadronic origin. Under the assumption that 
proton acceleration to high energies (~ 10 16 — 10 1 ' eV) 
is also viable in blazar jets, the production of charged pi- 
ons is a natural outcome of photopion (pn) interactions 
between the relativistic protons and the internally pro¬ 
duced synchrotron photons. The decay of charged pions 
results in the injection of secondary relativistic electron- 
positron pairs 2 . It is the synchrotron radiation of the 
latter that emerges in the GeV - TeV regime, in con¬ 
trast to the hadronic cascade scenario for blazar emission 
(Mannheim et al. 1991; Mannheim 1993), where the cas¬ 
cade emission mainly contributes to the 7 -ray regime 3 . 
As the synchrotron self-Compton emission from primary 
electrons may also emerge in the GeV - TeV energy band, 
the observed 7 -ray emission can be totally or partially 
explained by photohadronic processes, depending on the 
specifics of individual sources (Petropoulou et al. 2015). 
Since the luminosity of the pn component is directly con¬ 
nected to that of very high-energy (~ 2 — 20 PeV) neu¬ 
trinos, our approach allows us to associate the observed 
7 -ray blazar flux with the expected neutrino flux. The 
physical model we use has been described, in more gen¬ 
eral terms, by Dimitrakoudis et al. (2012) and Dimitrak- 
oudis, Petropoulou, & Mastichiadis (2014). 


3 MONTE CARLO SIMULATIONS 

In the first two papers of this series we presented the 
principles on which the BSV is built and mostly con¬ 
centrated on the statistical properties of blazars, such 
as distribution trends, average values of some important 
parameters, and compared those with observed distribu¬ 
tions and values in radio and X-ray surveys. In Paper 
III we tied our simulation to the absolute numbers of 
the Fermi- 2LAC catalogue and predicted the number of 

- We will commonly refer to them as electrons. 

3 The cascade is initiated by photon-photon absorption of 
very high-energy 7 -rays (7> TeV), which, in turn, are pro¬ 
duced by synchrotron radiation of secondary pairs from pion 
decays and Bethe-Heitler (pe) pair production. 


sources in the VHE band taking into account the extra- 
galactic background light (EBL) absorption. In Paper IV 
we calculated the integrated flux from the entire popula¬ 
tion of blazars of different types in the X-ray and 7 -ray 
bands. 

For the reader’s convenience we briefly summarise 
here our Monte Carlo simulations covering the radio 
through the 7 -ray bands, which are based on a number of 
ingredients, including: the (radio) blazar luminosity func¬ 
tion and evolution, a distribution of the Lorentz factor of 
the electrons and of the Doppler factor, a synchrotron 
model, an accretion disk component, the host galaxy, 
plus a series of 7 -ray constraints based on observed dis¬ 
tributions estimated using simultaneous multi-frequency 
data: the distribution of the ratio between high-energy 
and low-energy hump fluxes, the dependence of the 7 - 
ray spectral index on Op eSik , and that of the 7 -ray flux 
on radio flux density. Sources are classified as BL Lacs, 
FSRQ, or radio galaxies based on the optical spectrum, 
as in real surveys. Readers are referred to Paper I and 
II for full details. We stress that no assumption has been 
made in the simulations about the process responsible for 
the 7 -ray emission. 

As done in Paper III the SEDs were extrapolated 
to the VHE band by using our simulated Fermi fluxes 
and spectral indices and assuming a break at E = Eb rea k 
and a steepening of the photon spectrum by Ar. Our 
adopted values were Ebreak = 100 GeV and AP = 1 and 
Ebreak = 200 GeV and AP = 0.5, the latter being our 
default choice. VHE spectra were attenuated using recent 
estimates of the EBL absorption as a function of redsliift 
(Dominguez et al. 2011). 


3.1 The neutrino component 

The Monte Carlo simulations described above charac¬ 
terise fully the photon emission from the BL Lac popula¬ 
tion. To get to the neutrinos, we need to add a hadronic 
“prior”, which is built on the knowledge gained by fit¬ 
ting the SEDs of six individual BL Lac sources and their 
respective neutrino spectra in Petropoulou et al. (2015). 

We model the observed differential neutrino plus 
anti-neutrino (v + u) energy flux of all flavours (E„(E„)) 
as 

F V {E V ) = F 0 E- S exp J) , (1) 

where the parameters to be defined are Eo (char¬ 
acteristic energy), s (spectral slope), and Eo (normalisa¬ 
tion). 


3.1.1 Determination of Eo 

We set the characteristic energy Eo equal to the peak 
energy of the neutrino spectrum E„ iP , which, in good 
approximation, can be written as 
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Ev,p($, Z , i'peak) 
Eq — Ei/^p 


17.5 PeV / j _\ 2 /_ ^peak 
(I + 2) 2 1.10 j l 10 1S Hz 


( 2 ) 


where 8 is the Doppler factor and 2 is the source 
redshift (e.g. Dermer, Ramirez-Ruiz, & Le 2007). Here 
^peak is the observed, synchrotron peak frequency. This 
relation is valid as long as protons energetic enough to 
produce pions by interacting with the photons of the low- 
energy (synchrotron) hump of the SED can be accelerated 
in the jet. We assume that this condition is satisfied in 
all BL Lacs and discuss possible caveats in Sect. 5.2. 

Equation 2 requires as input 2 , <5, and iZp eak . These 
are all part of the output of the simulations so Eo can be 
easily derived source by source. 


3.1.2 Determination of s 

In Petropoulou et al. (2015) the neutrino spectra from six 
BL Lac objects were calculated self-consistently. It was 
shown that all neutrino spectra peaked at ~ E Vjp and 
E V F V {E U ) oc E~ s+1 with ( s ) ~ —0.35, namely the neu¬ 
trino spectra are relatively flat. This is not a matter of 
fine tuning or coincidence as, in our scenario, the number 
density of the target photons decreases with photon en¬ 
ergy 4 . Based on the range of values found by Petropoulou 
et al. (2015), the neutrino slope s was drawn from a Gaus¬ 
sian distribution centred at —0.35 with a = 0.12. 


3.1.3 Determination of Fo 
The normalisation Fq is given by 


F 0 


F u , tot E^p 1 


f°° dx x~ s e~ x ’ 


( 3 ) 


where x = E v / E VtP and x m j n is the minimum nor¬ 
malised neutrino energy. Because of the neutrino spec¬ 
tral shape, it is safe to set a: m m = 0 , and the integral 
reduces to the T function with argument 1 — s. In our ap¬ 
proach, the integrated neutrino flux is associated to the 
integrated 7 -ray flux as 


El,tot = Y V1 F 1 (> E~f ), (4) 

where we chose E 1 = 10 GeV 5 . By construction, 
Y vy encloses all the information on the optical depth 
for photopion interactions and the proton luminosity. 


4 This can be understood as follows: we assume that the pro¬ 
ton distribution extends up to 7 P ,max 7p,th; where 7 p ,th is 
the threshold energy for pir interactions with synchrotron pho¬ 
tons of frequency r'p eak . Protons with 7 P < 7 p ,thi which pion- 
produce on synchrotron photons with v > ^p eak , are respon¬ 
sible for the production of neutrinos with energies E v < E v ,p. 
Thus, the neutrino spectrum below its peak energy probes the 
soft part of the low-energy hump of the SED. 

5 Our results are totally independent of this choice. 


It has an upper limit, i.e. Y v - t ~ 3, which is obtained 
when synchrotron emission from p 7 r pairs accounts for 
the whole observed 7 -ray flux (e.g. Petropoulou & Mas- 
tichiadis 2015). To reach this result, we also made some 
simple assumptions about the pion production ratio and 
the energy of the leptons produced from the pion decay 
chain. However, there is no lower limit on Yl 7 : Yl 7 -C 1 
simply implies a leptonic origin for the 7 -ray blazar emis¬ 
sion. Application of the leptohadronic model to individ¬ 
ual BL Lacs with different properties, such as SED and 
redshift, resulted in Yl 7 values covering the range 0.1 — 2 
(Petropoulou et al. 2015). 

Based on the results of Petropoulou et al. (2015), and 
assuming that Yl 7 is the same independently of BL Lac 
class (but see Sect. 5.2.1) we adopt two different options 
for the ratio Vl 7 : 

(i) we assume a constant value for all BL Lacs, Yl 7 = 
0 . 8 , which is the average derived from the modelling of 
the six BL Lacs; 

(ii) we use a relation between Yl 7 and the observed 7 - 
ray luminosity above 10 GeV in erg s“ , Z/ 7 (> 10 GeV), 
i.e. logkl 7 = —0.58 x logL 7 (> 10 GeV) + 26.3, with 
Y v ., ^ 3. This is a simple fit to the values in Petropoulou 
et al. (2015) (see their Fig. 11) motivated by the possible 
trend found in the sample of six BL Lacs. In practice, 
this translates to Vl 7 = 3 for L 7 (> 10 GeV) < 3 x 10 44 
erg s _1 , that is for basically all LBL and ~ 80% of HBL. 

Combining all of the above, we derive the final ex¬ 
pression for the observed neutrino flux expected from 
each BL Lac object ( E 2 dN/dE units): 


EvF v (E„) = 

Y VJ FJ> 10 GeV) ( E„ \ ~ s+1 ( E„ \ (5) 

f x °°. dxx-e~* \E V J 6XP { Ev,p) 

J 'min 

The neutrino fluxes at various energies were then 
calculated based on this spectrum, which in turn depends 
on Ev, P , s, and Y„ 7 , and the 7 -ray flux and power from 
the simulations. Since Ev, p is fully determined by eq. 2 
and s covers a narrow range, we are left only with Y „ 7 as 
a possible “tuneable” parameter (see Sect. 5.5). 

Finally, the total NBG was computed as B = 
rFmax p av dp where ^ are the differential number 

dt di¬ 

counts and Fmi n and Fmax are the fluxes over which 
these extend. Since the IceCube data are “per neutrino 
flavour” our numbers, which refer to all neutrino flavours, 
were divided by three 6 . Ten simulations were run and the 
average was then calculated to smooth out the “noise” 
inherent to the Monte Carlo process. 


6 The neutrino flux produced at the source contains neutrinos 
of different flavours with an approximate ratio F Ve : F v ^ : 
F l ,_ =2:1:0. However, by the time they reach Earth their 
ratio will have changed to F v< , : E Uji : F lf _ =1:1:1 due 
to neutrino oscillations (Learned & Pakvasa 1995), as indeed 
observed (Aartsen et al. 2015). 
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Figure 1. The predicted neutrino (i/ + v) background per neu¬ 
trino flavour from BL Lacs. Different lines correspond to dif¬ 
ferent assumptions (starting from the top): Y„ 7 anti-correlated 
with L. y(> 10 GeV) for all BL Lacs (black dot short-dashed 
line) and HBL (black dot long-dashed line); Y|, 7 = 0.8 and 
®break = 200 GeV, AL = 0.5, for all BL Lacs (blue solid 
line) and HBL (blue dotted line); Y v ~f = 0.8 and -E brea k = 100 
GeV, AL = 1, for all BL Lacs (green long dashed line). All 
lines correspond to the mean value of ten different simula¬ 
tions. The (red) filled points are the data points from IceCube 
Collaboration (2014), while the open points are the 3(7 upper 
limits. The upper (magenta) short dashed lines represents the 
90% C.L. upper limits from The Pierre Auger Collaboration 
(2015) while the lower (cyan) short dashed line is the expected 
three year sensitivity curve for the Askaryan Radio Array (Ara 
Collaboration 2012). 


4 THE NEUTRINO BACKGROUND FROM 

BL LACS 

Fig. 1 shows our results, with different lines correspond¬ 
ing to different assumptions. Namely, and starting from 
the top: y „ 7 anti-correlated with L 7 (> 10 GeV) for all 
BL Lacs (black dot short-dashed line) and HBL (black 
dot long-dashed line); Y„ 7 = 0.8 and E brea k = 200 GeV, 
Ar = 0.5, for all BL Lacs (blue solid line) and HBL (blue 
dotted line); Y „ 7 = 0.8 and E bre ak = 100 GeV, AF = 1, 
for all BL Lacs (green long dashed line). The (red) filled 
circles are the data points from IceCube Collaboration 
(2014), while the open points are the 3cr upper limits, 
derived from the la ones by simply scaling them by a 
factor ~ 3.6 (Gehrels 1986). The upper (magenta) short 
dashed lines represents the 90% C.L. upper limits from 
The Pierre Auger Collaboration (2015) while the lower 
(cyan) short dashed line is the expected three year sensi¬ 
tivity curve for the Askaryan Radio Array (Ara Collabo¬ 
ration 2012 ). 

A few points can be made about this figure: 

(i) BL Lacs as a class can easily explain the whole 


NBC at high-energies (> 0.5 PeV) while they do not 
contribute much (~ 10%) at low-energies (< 0.5 PeV); 

(ii) HBL make up the bulk of the NBC in all cases 
up to ss 30 PeV, where they start to contribute less and 
less. This is where LBL take over, due to their lower i/p eak 
and therefore larger values of neutrino peak energy, given 
the inverse dependence between i/p eak and E„ iP (eq. 2 ). 
Note that HBL manage to dominate the neutrino out¬ 
put despite their small fraction (~ 5%) because of their 
relatively high 7 -ray, and therefore neutrino, fluxes and 
powers; 

(iii) there is very little difference between the _E brea k = 
200 GeV, AT = 0.5 and the E break = 100 GeV, Ar = 1 
cases. This is due to the fact that we relate the neutrino 
flux to E 7 (> 10 GeV) and most of the 7 -ray flux is at 
lower energies. Thus, the neutrino flux in the latter case 
is only slightly smaller than in the former; 

(iv) The difference between the constant and varying 
Y „ 7 case is small up to ~ 1 — 2 PeV. However, at 
E v > 5 PeV only the constant Y „ 7 case is fully consistent 
with the IceCube data, while the varying V „ 7 scenario 
is barely within the 3er upper limits. Moreover, current 
Auger upper limits already rule out our varying Y „ 7 sce¬ 
nario, while future ones will better probe the high-energy 
tail. The Askaryan Radio Array, by sampling the > 10 
PeV range with high sensitivity, will further constrain our 
model, together with the Antarctic Ross Ice Shelf An¬ 
tenna Neutrino Array (ARIANNA, Barwick et al. 2015), 
which will cover the 100 PeV ■ 10 EeV energy range; 

(v) if our model is correct, IceCube should see more 
E v > 1 PeV, and even E v > 2 PeV, events; alternatively, 
a non-detection will put strong constraints on Y ), 7 (see 
Sect. 5.5). 

Given the above discussion, from now on we refer to 
the Y „ 7 = 0.8 and E brea k = 200 GeV, AP = 0.5 case as 
our “benchmark” case. 

Giommi & Padovani (2015) have shown that if the 
Monte Carlo simulations are modified to include the 
strong dependence of v peak on radio power postulated by 
the blazar sequence (Fossati et al. 1998) the agreement 
with the observational data disappears, as the predicted 
7 -ray background above a few GeV turns out to be far in 
excess of the observed value. A similar thing happens for 
the predicted NBG, which, for example, turns out to be 
more than two orders of magnitude above the IceCube 
data at E v ~ 1 PeV. 

Fig. 1 gives a feeling for the average NBG but not for 
the dispersion intrinsic to our simulation process. This is 
shown in Fig. 2, where the results from the ten differ¬ 
ent simulations are displayed individually. Different lines 
correspond to different assumptions. Namely, our bench¬ 
mark case: average value (solid blue line) and individual 
runs (dotted blue lines); V „ 7 = 0.8 and E bre ak = 100 
GeV, Ar = 1: average value (long dashed green line) and 
individual runs (dot long-dashed green lines). The top 
line (cyan dot short-dashed line) represents the case with 
E brea k = 200 GeV, AF = 0.5 but Y v ~, = 3. As this is the 
maximum value expected in our model for this parameter 
(e.g. Petropoulou & Mastichiadis 2015) this curve repre¬ 
sents the largest NBG from BL Lacs. The black square is 
the neutrino flux of IceCube event 9 (Padovani & Resconi 
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Figure 2. The predicted neutrino (u+u) background per neu¬ 
trino flavour from BL Lacs for three different cases and show¬ 
ing the range of results from ten different simulations. Different 
lines correspond to different assumptions (starting from the 
top): y„ 7 = 3 (the maximum theoretical value: cyan dot short- 
dashed line) and -Ebreak = 200 GeV, AT = 0.5; Y v -y = 0.8 
and Ebreak = 200 GeV, AT = 0.5: average value (solid blue 
line) and individual runs (dotted blue lines); Y„.y = 0.8 and 
Ebreak = 100 GeV, AT = 1: average value (long dashed green 
line) and individual runs (dot long-dashed green lines). The 
(red) filled points are the data points from IceCube Collabo¬ 
ration (2014), while the open points are the 3cr upper limits. 
The black square is the neutrino flux of IceCube event 9, from 
Padovani & Resconi (2014), who have tentatively associated it 
with MKN 421, converted to “background” units (i.e. divided 
by 47r). 


2014) converted to “background” units (i.e. divided by 
47t). 

A few points can be made about this figure: 

(i) the dispersion amongst the ten different simulations 
for the two cases is typically not very large (< 0.3 dex), 
with only 1 — 2 of them being clear outliers. In all cases 
this is due to a single source, which has a neutrino flux 
~ 20 — 45 times higher than the source with the second 
largest flux, at the energy where a “bump” can be seen 
in the curves. This has to be expected given the nature of 
the simulations but these outliers happen only ~ 10—20% 
of the time; 

(ii) the presence of these outliers can explain at least 
some of the low-energy associations made by Padovani 
& Resconi (2014). The neutrino flux of IceCube event 9, 
which has been tentatively associated by these authors 
with MKN 421, is in fact fully compatible with one of 
the runs. Moreover, the number of events predicted by 
some of these outliers at low energies (E v < 2 PeV) is 
about twice as large as the mean of the ten simulations; 

(iii) the curve for the case with Y V1 = 3 defines the 
region of parameter space allowed by our model. At high 



Figure 3. The redshift distribution for the BL Lacs contribut¬ 
ing ~ 95% of the background associated with the benchmark 
case at 1 PeV (black long-dashed line) and for those with a 
measurable redshift (red solid line). The sources detectable by 
the 3FGL (blue dotted line) and 2FGL (magenta short-dashed 
line) catalogues are also indicated. 

energies this is well above the IceCube 3 a upper limits, 
which means that this upper bound is clearly ruled out 
by current data, and therefore that a leptonic component 
must be present in the y-ray emission of blazars. 


5 DISCUSSION 

5.1 Characterising the NBG sources 

As is normally the case for astrophysical backgrounds, 
most of the NBG is produced by a small fraction of the 
population. In our case, ~ 0.5% of the BL Lacs make 
~ 95% of the background associated with the benchmark 
case at 1 PeV, which is where our predictions are quite 
close to the IceCube data. We now focus on this sub¬ 
sample and describe some of its properties. 

Fig. 3 shows the redshift distributions for the whole 
sub-sample' (dashed line) and for the subset of BL Lacs 
with measurable redshift (red solid line). Note that the 
latter make up ~ 50% of the sub-sample, with the other 
50% consisting of BL Lacs with no measurable redshift. 
These are BL Lacs whose maximum equivalent width 
(EW) is < 2 A, or for which the non-thermal light is 
at least a factor 10 larger than that of the host galaxy, 
and are therefore deemed to have a redshift which can¬ 
not be typically measured (Paper I). As expected, the two 
distributions are quite different with (z) ~ 1.3 and ~ 0.9 
respectively. This is due to the fact that the sources with 

7 We note that this is not significantly different from the red¬ 
shift distribution of the entire sample. 
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no measurable redshift are mostly those with the optical 
spectra swamped by non-thermal emission, which tend 
to have a more powerful jet and therefore are on average 
at larger intrinsic redshift. 

The fraction of sources in the sub-sample detectable 
by the 3FGL Fermi catalogue (The Fermi-LAT Collab¬ 
oration 2015) is ~ 11%. However, these make ~ 50% of 
the overall background associated with the benchmark 
case, i.e. ~ 15% of the NBC. We also find that ~ 7% of 
the sources would be detectable by the 2FGL (Nolan et 
al. 2012), making ~ 40% of the background associated 
with the benchmark case, i.e. ~ 12% of the NBG. These 
numbers are consistent with the preliminary results of 
Gliisenkamp et al. (2015), who find no evidence of neu¬ 
trino emission and a maximal contribution from Fermi 
2LAC (Ackermann et al. 2011) blazars ~ 20% 8 . 

The (intrinsic) redshift distributions of the sources 
detectable by the 3FGL (blue dotted line) and 2FGL 
(magenta short-dashed line) catalogues are also shown 
in Fig. 3. 

5.2 Possible caveats 

5.2.1 LBL 

The inclusion of the hadronic “prior” into the Monte 
Carlo simulations is based upon the results derived from 
individual SED fitting of HBL (Pctropoulou et al. 2015). 
The application of our model to LBL, therefore, could 
not be straightforward. Thus, we discuss here possible 
caveats. 

A direct application of eq. 2 to LBL requires that: 
(1) the photons of the low energy SED component are the 
main targets for photopion interactions; ( 2 ) acceleration 
of protons to ultra-high energies (UHE), e.g. 10 19 — 10 2 ° 
eV, takes place in the jets of these sources. These assump¬ 
tions, in turn, suggest that specific conditions should 
prevail in the emitting region of individual LBL, in or¬ 
der to explain their 7 -ray emission in terms of photo- 
hadronic processes. In particular, if the 7 - 1 'ay emission 
is synchrotron radiation of UHE secondary pairs (simi¬ 
larly to HBL), then the emitting region should contain 
very weak magnetic fields, e.g. B <C 3/rG (see also eq. 9 
in Petropoulou & Mastichiadis 2015). Alternatively, the 
7 -ray emission of LBL could be the result of a hadronic 
cascade (Mannheim et al. 1991; Mannheim & Biermann 
1992). If the above conditions can be realised in the jets 
of LBL, then application of eq. 2 shows that the LBL 
contribution to the NBG takes place at high energies, i.e. 
> 100 PeV (see Fig. 1). 

One could relax, however, assumption (2), if the tar¬ 
gets for the photopion interactions were synchrotron pho¬ 
tons with frequencies above = 10 13 Hz (see generic 
SED in Fig. 8 in Petropoulou & Mastichiadis 2015). In 
this case, application of eq. 2 to LBL would also result 
in the production of neutrinos with energies similar to 
those expected from HBL 9 . In this scenario, the neutrino 

8 This is however derived assuming an V/ 2-46 spectrum, 
which is much softer than our predictions. 

9 A preliminary application to Ap Librae results in ~ 0.5 

and a neutrino peak energy of several PeV. 


emission from both HBL and LBL sources is expected to 
have a similar peak energy. 

To test the robustness of the results shown in Figs. 1 
and 2 we therefore artificially shifted the iZpeak distribu¬ 
tion of LBL so that it had the same mean as HBL, using 
the same value of Y vy . We then repeated the Monte Carlo 
simulations and found that the total NBG flux increased 
only by ~ 15% overall. This is due to the fact that LBL 
have much smaller 7 -ray, and therefore neutrino, fluxes 
than HBL. It then follows that our results on the NBG 
do not depend on a detailed modelling of LBL and are 
therefore robust. 


5.2.2 Energetics 

Leptohadronic models are known to predict relatively 
high jet powers, driven by their relativistic proton content 
(e.g. Zdziarski & Bottcher 2015; Petropoulou et al. 2015, 
and references therein). Indeed, the estimated jet power 
for the six HBL to which Petropoulou et al. (2015) ap¬ 
plied their leptohadronic model, and which lie at the basis 
of our calculations, is quite high (Lj et ~ 2 x 10 48 — 8 x 10 49 
erg s ' 1 [see their Tab. 4]). 

One could compare L ]et with the Eddington lumi¬ 
nosity, LEdd, of those sources. Assuming a typical black 
hole mass of 3 x 10 s Mg (Plotkin et al. 2011), this would 
be ~ 3.8 x 10 46 erg s _1 , which, in turn, would imply that 
the jet power is > 50 (and up to ~ 2, 000) times larger. 
These large Ljet/L/Edd ratios are not necessarily alarm¬ 
ing, for the following reasons: 1. out of the six BL Lacs 
modelled in Petropoulou et al. (2015) only for MKN 421 
is there an estimate for the black hole mass (Sbarrato et 
al. 2012 ) (which happens to be equal to the value we have 
assumed), which means that a comparison against an un¬ 
certain LEdd is not very informative; 2. Ghisellini et al. 
(2014), by applying a leptonic model to a large sample of 
mostly FSRQ, have derived at low powers Lj e t ~ I/Edd, 
with some sources reaching Ljet/L/Edd ~ 30. At least some 
of the sources studied by Petropoulou et al. (2015) (in¬ 
cluding MKN 421) lie close to the upper range of that 
distribution; 3. the total jet power does not need to be 
limited by LEdd, as is observed, for example, in the case 
of gamma-ray bursts, which are super-Eddington by huge 
amounts (« 10 10 ; e.g. Piran 2004). 

A more meaningful comparison would be that of Lj e t 
with Me 2 , where M is the accretion rate. One in fact ex¬ 
pects Ljet ~ tjMc 2 with ej < 1.5 (Zdziarski & Bottcher 
2015, and references therein). In the case of radiatively 
efficient accretors, the accretion rate is normally derived 
from the relation Me 2 = Ldisk/??, where 7 is the radiative 
efficiency. Ldisk, in turn, is estimated through the broad 
emission lines, which are used, via a template, to derive 
the luminosity of the entire broad line region (Lblr). The 
latter is therefore a proxy for the accretion disk luminos¬ 
ity, modulo a covering factor, which is usually taken to 
be ~ 0.1 (which means Ldisk ~ IOLblr: Ghisellini et 
al. 2014, and references therein). None of the six sources 
studied by Petropoulou et al. (2015) shows broad lines, 
while only two sources display narrow lines in their spec¬ 
tra. In fact, our sources are most likely LERGs and there¬ 
fore radiatively inefficient accretors. In short, it is not 
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possible to derive the accretion power of these sources, 
since they do not exhibit broad line features. Most impor¬ 
tantly though, we could not relate Ldisk to Me 2 because 
the simple relation, which is valid for radiatively efficient 
accretors, breaks down for LERGs (Sbarrato, Padovani, 
& Ghisellini 2014, and references therein). The point re¬ 
mains that we do not expect high values of M, as other¬ 
wise these sources could not be LERGs. 

If our model will be confirmed by future IceCube 
detections, this will simply mean that the current picture 
of accretion in broad-lined, efficient accretor AGN cannot 
be applied to HBL. 

As a final point, we note that the energetics for BL 
Lacs in the proton synchrotron scenario are more rea¬ 
sonable than in the leptohadronic one (Dimitrakoudis, 
Petropoulou, & Mastichiadis 2014) but in that case the 
neutrino flux from BL Lacs would peak at E u ~ 200 PeV 
and would be much lower than in our model. 

It could be argued that, given the high powers in rel¬ 
ativistic protons, proton-proton (pp) collisions with the 
background (thermal) protons in the inner jets of blazars 
could result in a substantial production of neutrinos. 
While a detailed study of this process goes beyond the 
scope of this paper, we note that, having a neutrino flux 
from the pp mechanism equal to that from the pn interac¬ 
tion at ~ 1 PeV, would require thermal proton number 
densities « 10 8 cm -3 . These, in turn, would translate 
into a (thermal) proton energy density ~ 1 — 2 orders of 
magnitude larger than that derived by Petropoulou et al. 
(2015) and a correspondingly larger Ljet- We thus con¬ 
tend that this scenario is not plausible (see also Atoyan 
& Dermer 2003). 


5.2.3 FSRQ 

The model we used for including neutrino emission in the 
Monte Carlo simulation has not been applied to FSRQ, 
which we therefore excluded from our calculations. A 
proper treatment of the neutrino emission from individ¬ 
ual FSRQ is required in order to make robust predictions 
about their cumulative contribution, but this lies outside 
the scope of the present study. Recently, Dermer, Murase, 
& Inoue (2014) and Murase, Inoue & Dermer (2014) have 
calculated the neutrino background from FSRQ based on 
the blazar sequence and assuming a leptonic origin of the 
blazar 7 -ray emission. They have taken into account pho- 
topion interactions with the non-thermal emission from 
the jet and the external photon fields (BLR and dusty 
torus). The latter have been shown to provide most of 
the contribution to the total neutrino output. 

Here instead, we calculated the NBG from FSRQ 
assuming that they fall within the same scenario of BL 
Lacs, an assumption that has no astrophysical basis. We 
find that their contribution to the NBG is somewhat en¬ 
ergy dependent but overall only « 30% that of BL Lacs. 
We caution the reader that this estimate neglects the role 
of external photon fields in the photopion interactions 
and, therefore, in the final neutrino output. As such, this 
contribution should be considered as a lower limit. 


5.2.4 Masquerading BL Lacs 

Paper I and II discussed the existence of sources, which 
appeared “BL Lac-like” only because their emission lines 
were heavily diluted by strong non-thermal emission. 
These objects, which are therefore “masquerading” BL 
Lacs and are in fact misclassified FSRQ, make up a very 
small fraction (~ 2%) of our BL Lacs. However, as is the 
case for HBL, because of their relatively high 7 -ray, and 
therefore, neutrino fluxes, their contribution to the NBG 
is « 60% that of all BL Lacs. Our calculations on the neu¬ 
trino emission were made under the assumption that only 
photons produced internally in the jet, i.e. synchrotron 
photons, are the targets for photopion interactions. In 
masquerading BL Lacs any line- or blackbody-likc emis¬ 
sion external to the jet is hidden below the non-thermal 
synchrotron emission that is produced within the jet. We 
argue that our model is also applicable to these sources, 
as long as u' syn > it' x , where u' syn and u' x ~ r 2 u e x (where 
r is the Lorentz factor) are the energy densities of syn¬ 
chrotron and external photons, respectively, as measured 
in the rest-frame of the emitting region. Obviously, our 
argument becomes even stronger if the boosting of the 
external energy density becomes less efficient (e.g. the 
non-thermal emitting region is located much further out 
than the external photon field region). 


5.2.5 Flaring sources 

One of the defining characteristics of blazars is their large 
variability known to occur in all parts of the electromag¬ 
netic spectrum. The specific amount depends on the en¬ 
ergy band and on the blazar type, with values ranging 
from a factor of a few in the radio band up to a fac¬ 
tor 10,000 at GeV energies (e.g. Aharonian et al. 2007; 
Giommi 2015). The amount of variability of the neu¬ 
trino flux in blazars and how this correlates with 7 -ray 
variability is not known. As the number of astrophysi¬ 
cal neutrino events detected so far is still very small, one 
could speculate that large flaring events from one or a few 
of the brighter sources in Fig. 2 may produce a fair frac¬ 
tion of the observed events. This is not straightforward 
though, and requires detailed modelling of flaring sources 
(e.g. Reimer et al. 2005). A more recent application to the 
flaring MKN 421 also suggests that an accumulation of 
flaring events from an individual source may be needed in 
order to have a statistically significant neutrino detection 
(Petropoulou et al., in preparation). Here we consider the 
emission arising from the entire population of BL Lacs 
and assume that flares from single sources will be diluted 
by the integrated emission from all blazars. 


5.3 Comparison with previous results 

The idea of blazars being sources of high-energy neutri¬ 
nos dates back to long before the detection of sub-PeV 
neutrinos with IceCube (IceCube Collaboration 2013) 
and has since been explored in a number of studies (e.g. 
Mannheim 1995; Halzen & Zas 1997; Miicke et al. 2003; 
Tavecchio & Ghisellini 2015; Kistler, Stanev, & Yiiksel 
2014; Dermer, Murase, & Inoue 2014; Murase, Inoue & 
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Dermer 2014). Our approach has some similarities but 
many differences with previous work, as detailed below. 

5.3.1 Similarities 

(i) the BL Lac 7 -ray emission has a (photo)hadronic 
origin (at least for the models presented in Fig. 4); 

(ii) in BL Lacs the targets for photopion interactions 
are the low-energy synchrotron photons. 

5.3.2 Differences 

(i) we use as a starting point the knowledge gained 
from detailed SED fitting of BL Lacs instead of using a 
generic neutrino spectrum (e.g. Mannheim, Protheroe & 
Rachen 2001; Kistler, Stanev, & Yiiksel 2014; Murase, 
Inoue & Dermer 2014). By establishing a connection be¬ 
tween the 7 -ray and neutrino emission for each source 
(see eqs. 2-5), we are able to assign to each simulated 
BL Lac in the Monte Carlo code a unique neutrino spec¬ 
trum. We then calculate the NBG by summing up the 
fluxes of all sources in each energy bin; 

(ii) for the calculation of the NBG we do not normalise 
a priori a generic neutrino spectrum to the extragalactic 
7 -ray background (EGB) (e.g. Mannheim 1995; Miicke et 
al. 2003). In fact, we do not need to, as our simulation 
naturally reproduces the observed EGB above 10 GeV 
(Giommi & Padovani 2015); 

(iii) the NBG spectrum is not a priori normalised to 
the IceCube observations (e.g. Tavecchio & Ghisellini 
2015). Instead, for a specific choice of Y v - n which is the 
only tuneable parameter in our framework, we compare 
our model predictions with the IceCube data; 

(iv) the maximum proton energy is taken to be a few 
times larger than the threshold energy for photopion in¬ 
teractions with the peak energy synchrotron photons of 
the low-energy hump. This is usually lower than the val¬ 
ues used in previous studies (e.g. Halzen & Zas 1997; 
Miicke et al. 2003), which also explains the difference in 
the peak energies of the NBG; 

(v) the 7 -ray emission of individual BL Lacs in our ap¬ 
proach is a combination of synchrotron radiation emitted 
by electron-positron pairs produced through ^ decay 
and synchrotron self-Compton from primary electrons. 
The cascade emission initiated by 7 r° 7 -rays has a negli¬ 
gible effect in the formation of the blazar SED. This is 
in contrast to previous studies, where the blazar 7 -ray 
emission is explained either as proton synchrotron radia¬ 
tion (e.g. Miicke et al. 2003) or as cascade emission (e.g. 
Halzen & Zas 1997; Kistler, Stanev, & Yiiksel 2014). 


5.3.3 Detailed comparison 

Fig. 4 compares the predicted neutrino background for 
our benchmark case for all BL Lacs (blue solid line) and 
HBL (blue dotted line) with some of the previous results. 
I 11 chronological order, these are: Mannheim (1995) (long 
dashed cyan line, upper limits at low energies), Halzen & 
Zas (1997) (short dashed green line), Miicke et al. (2003) 
(dot long-dashed black lines), and Tavecchio & Ghisellini 
(2015) (dot short-dashed magenta line). The two curves 



Figure 4. The predicted neutrino background per neutrino 
flavour for Y „. Y = 0.8 and Fbreak = 200 GeV, AT = 0.5, 
for all BL Lacs (blue solid line) and HBL (blue dotted line) 
compared to previous results. Namely, in chronological or¬ 
der: Mannheim (1995) (long dashed cyan line; upper limits at 
low energies), Halzen & Zas (1997) (short dashed green line), 
Miicke et al. (2003) (dot long-dashed black lines: LBL, upper 
curve; HBL, lower curve), and Tavecchio & Ghisellini (2015) 
(dot short-dashed magenta line). The (red) filled points are 
the data points from IceCube Collaboration (2014), while the 
open points are the 3<r upper limits. See text for details. 

from Miicke et al. (2003) represent the maximum contri¬ 
bution expected from LBL (upper curve) and HBL (lower 
curve), respectively. A few things about Fig. 4 are worth 
mentioning: 

(i) the model by Mannheim (1995) at first glance is the 
one that best describes the IceCube data. This, taken at 
face value, would imply that radio-loud AGN explain the 
entire NBG, something that contradicts the preliminary 
IceCube results of Gliisenkamp et al. (2015), who find 
a maximal contribution from Fermi 2LAC (Ackermann 
et al. 2011) blazars ~ 20%. However, since it gives only 
upper limits at low energies, it could be still reconciled 
with the data. This model has a very different shape as 
compared to the others because it includes two hadronic 
components, i.e. a low-energy soft one (E u < 2 PeV), pro¬ 
duced through pp collisions of the escaping CRs from the 
blazar jet with the ambient medium, and a high-energy 
flat one (E v > 2 PeV), related to pn interactions of CRs 
with the synchrotron photons in the blazar jet; 

(ii) the model by Halzen & Zas (1997), although very 
close to ours at low energies, lies above the 3 a upper 
limits at E v > 5 PeV, while the sum of the two curves 
by Miicke et al. (2003) remains consistently below the 
IceCube data. Although the model curve of Tavecchio & 
Ghisellini (2015) passes through the data points, this is 
by construction, i.e. the NBG was a priori normalised to 
the IceCube data. Moreover, this model might also con- 
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tradict the IceCube results of Gliisenkamp et al. (2015) 
mentioned above (pending the different spectral shapes); 

(iii) apart from the Tavecchio & Ghisellini (2015) 
model and the low-energy (E„ < 2 PeV) part of the 
Mannheim (1995) model, the NBG spectrum below its 
peak is relatively hard in all other models shown in Fig. 4. 
This is related to the assumed proton power-law index 
and the spectrum of target photons, which are both sim¬ 
ilar among the models. The diversity of the peak energy 
of the NBG spectrum, on the other hand, reflects mainly 
the different model assumptions on the maximum proton 
energy in blazars. 


5.4 Cosmic rays 


Neutrons that are produced in photopion interactions 
are an effective means of CR escape from the system. 
They are unaffected by its magnetic field, their decay 
time is high enough to allow them to escape freely be¬ 
fore reverting into protons (Kirk & Mastichiadis 1989; 
Begelman, Rudak, & Sikora 1990; Giovanoni & Kazanas 
1990; Atoyan & Dermer 2003), and they are unaffected by 
adiabatic energy losses that the protons may sustain in 
the system as it expands (Rachen & Meszaros 1998). For 
the purposes of this paper, we will assume that “neutron 
conversion” is the injection mechanism of CRs from BL 
Lacs into the intergalactic medium (e.g. Kistler, Stanev, 
& Yiiksel 2014), while we neglect any additional contri¬ 
bution to the neutrino and CR fluxes from direct proton 
escape (e.g. Essey et al. 2010; Kalashev, Kusenko, & Es- 
sey 2013). This is a valid assumption as long as the es¬ 
caping protons are susceptible to adiabatic energy losses 
(e.g. the emitting region lies within an expanding jet) 
(Rachen & Meszaros 1998) and may end up carrying a 
negligible fraction of the UHECR flux. In this regard, the 
CR flux estimates that follow should be considered as a 
lower limit. 

The neutron spectrum, and thus, the escaping CR 
spectrum from each BL Lac, can be directly related to 
the neutrino spectrum (see e.g. Fig. 9 in Dimitrakoudis 
et al. 2012), at least in the optically thin regime for 
pn interactions 10 . Under certain simplifying assumptions, 
namely: (1) E n ~ 20 E v ; (2) E p ~ E n - and (3) pro¬ 
duction of one 7v± pair per interaction, which leads to 
E n dN/dE n « (l/6)E l ,dN/dE l , (see also Kistler, Stanev, 
& Yiiksel 2014), we may write 


E: 


. dN 
' dE p 


20 2 dN 
6 v dE v 


( 6 ) 


Having already calculated the NBG flux in Sect. 4, 
we can easily apply the above relation to all BL Lacs and 
estimate their contribution to the injected CR spectrum 
without taking into account any propagation effects. This 
will be a “copy” of the total NBG spectrum shown in 
Figs. 1 and 2, translated to higher energies by a factor of 
~ 20 and scaled in flux by a factor 11 of 10. 


10 The BL Lac emitting region is optically thin to p7r interac¬ 
tions with typical optical depths r p7r ~ 10 -6 — 10 —4 (see e.g. 
Atoyan & Dermer 2001; Petropoulou et al. 2015). 

11 Equation 6 was derived taking into account neutrinos of 


The resulting CR spectrum from HBL for our bench¬ 
mark case peaks at E p ~ 2.5 x 10 18 eV with a flux 
~ 3.0 x 10 24 eV 2 m -2 s” 1 sr _1 , only slightly above the 
observational data (see High Resolution Fly’S Eye Col¬ 
laboration et al. 2009; Abu-Zayyad et al. 2013; The Pierre 
Auger Collaboration 2013). A proper comparison would 
require the calculation of the expected proton fluxes at 
Earth taking into account propagation effects. Namely: 
(i) photopion (e.g. Mticke et al. 1999) and photopair 
(Blumenthal 1970) production on the background photon 
fields (cosmic microwave, cosmic infrared, and cosmic op¬ 
tical backgrounds); (ii) adiabatic energy losses; and (iii) 
magnetic deflections due to the intergalactic magnetic 
field. 

From the above, it becomes evident that a detailed 
treatment of LBL is also required, since the specifics of 
their SED modelling (e.g. the maximum proton energy) 
will affect our predictions of the UHECR spectrum. We 
plan to calculate the CR flux predicted by our scenario 
in detail in a future paper. 


5.5 Model predictions 

Our model makes specific predictions on the detectability 
of E„ > 2 PeV events, that is above the current maximum 
energy of IceCube neutrinos. We have used the effective 
areas from IceCube Collaboration (2013) and our NBG 
to predict the number of 2 < E v < 10 PeV neutrino 
events IceCube should see (as 10 PeV is the maximum 
energy for which IceCube effective areas are available). 
We get AC ~ 4.6 events for our benchmark case assum¬ 
ing the Glashow resonance (Glashow 1960) is not relevant 
since we are dealing with proton-photon interaction (An- 
chordoqui et al. 2005) (and AC ~ 7 otherwise). For our 
other case (Y„ 7 = 0.8 and E brea k = 100 GeV, AF = 1.0) 
we derive AC ~ 4.0 (or 6.1) events. Since the NBG of 
our models peak at E v > 10 PeV by making an educated 
guess on the effective areas we estimate an additional 2—3 
events up to ~ 100 PeV. Noting that the 3<r upper limit 
for 0 events is 6.6 (Gehrels 1986), we are actually close to 
being inconsistent with the IceCube non-detections. How¬ 
ever, our number is probably biased because Y„ 7 = 0.8 
is likely an upper limit. This was derived, in fact, from a 
small sample of BL Lacs, which might represent the tip 
of the iceberg in terms of neutrino emission, as they were 
selected as the most probable candidates. For example, 
if Y„ 7 = 0.3 then for our benchmark case AC « 3 (4) 
for 2 < E v < 100 PeV, which is well within the 2 a limit 
for 0 events. In any case, our predictions should easily be 
testable by IceCube in the next few years. 

Turning things around, we showed that the shape 
and the peak energy of the neutrino spectrum expected 
from a single BL Lac are mostly determined by its low- 
energy emission and the Doppler factor of the source (see 
Sect. 3.1). Thus, the calculation of the diffuse neutrino 
emission from BL Lacs can be considered to be robust. 
Most importantly though, our method has only one tune¬ 
able parameter, namely the ratio Y„ 7 (eq. 5). A com¬ 


all-flavours, whereas the NBG flux shown in Figs. 1 and 2 is 
per neutrino flavour. 
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Figure 5. The electromagnetic and neutrino extragalactic 
backgrounds predicted by our simulations in the energy range 
100 MeV - 300 PeV. The left side of the plot (E < 1 TeV) 
shows the 7 -ray background for various classes compared to 
the total extragalactic electromagnetic emission observed by 
Fermi -LAT (adapted from Paper IV), whereas the right side 
(E > 10 TeV) illustrates our prediction for the neutrino back¬ 
ground (all flavours) for our benchmark case (-Ebreak = 200 
GeV, AT = 0.5) for all BL Lacs (blue solid line) and HBL (blue 
dotted line) and Y v 7 ranging between 0.8 (upper curves) and 
0.3 (lower curves; see text for details). The (red) filled points 
are the (all flavours) data points from IceCube Collaboration 
(2014), while the open points are the 3<r upper limits. 

parison of the model predicted NBG with current Ice- 
Cube upper limits and, ultimately, future detections at 
Eu > 2 PeV, can be used to constrain the value of y„ 7 . 
In other words, this would provide an indirect way of 
probing the origin of the BL Lac 7 -ray emission. 

5.6 The big picture 

Fig. 5 displays both the electromagnetic and neutrino 
extragalactic backgrounds predicted by our simulations 
and the available measurements in the energy range 100 
MeV - 300 PeV. The left side shows the 7 -ray background 
compared to the total extragalactic electromagnetic emis¬ 
sion observed by Fermi -LAT (adapted from Paper IV), 
whereas the right side illustrates our prediction for the 
NBG (all flavours) for our benchmark case for all BL 
Lacs (blue solid line) and HBL (blue dotted line) and 
y „ 7 ranging between 0.8 (upper curves) and 0.3 (lower 
curves), the latter value being more consistent with the 
IceCube high-energy non-detections. 

The EGB can be approximated by a power law with 
exponential cutoff having Y ~ 2.3 and a break energy 
~ 280 GeV (Ackermann et al. 2015), the latter very likely 
due to the EBL absorption of 7 -ray photons from dis¬ 
tant (z > 0.3) sources (e.g. Ajello et al. 2015, and Paper 
IV). As a simple extrapolation of the EGB power law to 
the PeV energy range goes through the IceCube data, it 
might be tempting to assume that there is a single class 


of sources that explains both the EGB at E < 10 GeV 
and the NBG below ~ 0.5 PeV. This population cannot 
be the blazar one, for the following two reasons: (i) in 
the BSV scenario, blazars contribute ~ 50% — 70% to 
the total EGB at E < 10 GeV, while BL Lacs may ex¬ 
plain almost 100% of the EGB flux at E > 100 GeV; (ii) 
similarly, BL Lacs contribute only ~ 10% to the NBG 
at energies < 0.5 PeV, while they may fully explain the 
observed NBG above 0.5 PeV. If starburst galaxies can 
explain part of the EGB at E > 100 MeV (e.g. Lacki, Ho- 
riuchi, & Beacom 2014, and references therein) then they 
could also be a promising candidate class for explaining 
the sub-PeV IceCube neutrinos (e.g. Loeb & Waxman 
2006; Stecker 2007). We note, in fact, that in proton- 
proton scenarios of 7 -ray emission, relevant to starburst 
galaxies, the neutrino and 7 -ray spectra have the same 
power law index as the parent proton population (e.g. 
Kelner, Aharonian, & Bugayov 2006). 

Alternatively (or at the same time) the low-energy 
neutrino events could also have a Galactic component 
(e.g. Padovani & Resconi 2014). In any case, if there is 
a different class of sources contributing to the sub-PeV 
energy range, there is still room for individual BL Lac 
sources, like MKN 421 (see Fig. 2 and relevant discus¬ 
sion). Finally, we note that the EGB in Fig. 5 shows the 
sum of unresolved and resolved 7 -ray emission of the ex¬ 
tragalactic sky, whereas in the case of IceCube neutrinos 
we are not yet in the position to distinguish between a 
resolved and an unresolved contribution. The current sta¬ 
tus of neutrino astronomy, therefore, somewhat resembles 
that of 7 -ray astronomy in its very early days (i.e. those 
of SAS-2 and COS-B). 

The scenario, which appears to emerge by comparing 
our model NBG with the data is the following: at low en¬ 
ergy (Eu < 0.5 PeV) BL Lacs can only explain ~ 10% of 
the IceCube data. Some other population/physical mech¬ 
anism needs to provide the bulk of the neutrinos. How¬ 
ever, this does not exclude the possibility that individual 
BL Lacs still make a contribution at the « 20% level 
to the IceCube events. At high energy (Eu > 0.5 PeV) 
BL Lacs can account fully for the IceCube data. The 
strong implications of our scenario are: 1. IceCube should 
soon start resolving at least some of the NBG; 2. IceCube 
should also detect events at Eu > 2 PeV in the next few 
years. 


6 CONCLUSIONS 

We have included in the blazar simplified view scenario, 
which reproduces extremely well the statistical properties 
of blazars from the radio to the 7 -ray band, a hadronic 
component and calculated via a leptohadronic model the 
neutrino background produced by the whole BL Lac 
class. For the first time, this is done by summing up the 
fluxes of all the BL Lacs generated by a Monte Carlo sim¬ 
ulation, each with its own different properties. Our main 
results can be summarised as follows: 

(1) BL Lacs as a class can easily explain the whole 
neutrino background at high-energies (> 0.5 PeV) while 
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they do not contribute much (~ 10 %) at low-energies 
(< 0.5 PcV). 

(2) Individual BL Lacs, however, including some of the 
sources selected as possibly associated with some IceCube 
events by Padovani & Resconi (2014), can still make a 
contribution at the « 20 % level to the low-energy events. 

(3) Given our Monte Carlo approach, we can easily 
characterise the BL Lacs producing most of the neutrino 
background. These are of the HBL type (up to E v « 
30 PeV), have relatively low redshifts ((z) ~ 1-3), and 
about half of them have their optical spectra swamped 
by non-thermal emission and therefore an unmeasurable 
redshift. The (small) fraction of Fermi detectable sources 
in our scenario is consistent with a preliminary IceCube 
likelihood analysis on Fermi blazars. 

(4) A simultaneous look at the 7 -ray and neutrino 
backgrounds leads us to suggest that another popula¬ 
tion/physical mechanism could explain both the former 
at E < 10 GeV (since blazars dominate above that en¬ 
ergy) and the latter at E„ < 0.5 GeV. This might include 
star forming galaxies, although a Galactic component for 
the low-energy IceCube events could also be possible. 

(5) Our simulations at face value predict that IceCube 
should soon start resolving a fraction of the neutrino 
background and detect events at E u > 2 PeV. A non¬ 
detection will constrain the value of our only tuneable 
parameter (Y„ 7 ) and will provide an indirect way of prob¬ 
ing the origin of 7 -ray emission in BL Lacs. 
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